function   [simulation_results, h] = generate_recovery_dynamics(S, eta_adjustment, N_output )
% h is a handle to the figure of the post-crisis GDP response to credit spread
% simulation_results is a struct of simulation results.
% If N_output==1, then we only want to get the post-crisis GDP response and
% simulation results. Otherwise other output will be printed to pdfs.
mainpath = '/Users/wenhao/Dropbox/Research/Systemic Risk Project/Model_Only w dimension' ;
cd(mainpath);
addpath(genpath( 'Codes_Common Components'))   % Add other commonly shared files to the path. 

if( nargin==0 )
    S = load('solutions/baseline_rational.mat');  eta_adjustment = 0.05; N_output=1;
elseif( nargin==1 )   % By default only output the first figure
    eta_adjustment = 0.05;  N_output=1;
elseif( nargin==2 )
    N_output=1;  
end

%% Part 0: Load solutions
if(N_output==1)
    T_sim = 10000;
else
    T_sim = 50000;
end
simulation_results = simulation(  S , eta_adjustment, T_sim   );

N_sim = simulation_results.N_sim;
dt = simulation_results.dt;
dN_sim = simulation_results.dN_sim;
w_sim = simulation_results.w_sim;
w_sim_ode45 = simulation_results.w_sim_ode45;
productivity_sim = simulation_results.productivity_sim ;
credit_spread_sim = simulation_results.credit_spread_sim ;
xK_sim = simulation_results.xK_sim;
credit_sim = xK_sim.*w_sim;

CrisisStart = simulation_results.CrisisStart;
InCrisis = simulation_results.InCrisis;
RecessionStart = simulation_results.RecessionStart;
InRecession = simulation_results.InRecession;

%% Part 1: Match GDP impulse resonse by local projection method
% Calculate the Impulse Response of GDP to credit spread. 
lag_year_max = 5;
lag_years = (0:0.25:lag_year_max) ;   % quarterly
InBuffer = ( (1:N_sim)>(lag_year_max/dt+1)  &  (1:N_sim)<(N_sim-lag_year_max/dt-1) )';   % Make sure my calculation is not out of the data limit
types = ["Normal",  "Crisis", "Recession"];
coef_pre_all = cell(3,1);
coef_post_all =  cell(3,1);
for( k = 1:numel(types) )
    type =types(k);
    if( type=='Crisis' )    
        Start_Index = find(  CrisisStart  &  InBuffer  );
    elseif( type=='Recession' )
        Start_Index = find(  RecessionStart  & ~CrisisStart  &  InBuffer );
    elseif( type=='Normal' )
        Start_Index = find( ~RecessionStart & ~CrisisStart & InBuffer & (1:N_sim)'< 10^5  );   % Restrict the amount of calculation.
    end
    pre_spread = credit_spread_sim( Start_Index-1 );
    post_spread =  credit_spread_sim( Start_Index );
    coefs_pre = zeros(  numel(lag_years), 1 );  % Coefficients of regressing GDP to pre-crisis spread
    coefs_post = zeros(  numel(lag_years), 1 ); % Coefficients of regressing GDP to post-crisis spread
    for( i=1:numel(lag_years) )  % Lag of credit spread and GDP
        lag_year=lag_years(i);
        GDP_lag =  productivity_sim( Start_Index+lag_year/dt );
        b= [  ones( numel(GDP_lag), 1 )   pre_spread   post_spread ] \ GDP_lag  ;
        coefs_pre(i) = b(2) ;    
        %b= [  ones( numel(GDP_lag), 1 )      post_spread ] \ GDP_lag  ;
        coefs_post(i) = b(3) ;
    end 

    coef_pre_all{k} = coefs_pre;
    coef_post_all{k} = coefs_post;
end

% Compare the Impulse Responses of GDP to credit spread in one figure for all cases. 
types = ["Pre",  "Post"];
for( k=1:numel(types) )
    type = types(k);
    if( type=='Pre' )    
        coef =  coef_pre_all;
    else
        coef =  coef_post_all;
    end
    q = sum(  credit_spread_sim<mean( credit_spread_sim ) ) / numel( credit_spread_sim );
    change = 3.5*std(credit_spread_sim) ;   % 3.5 sigma of credit spread.
    ratio = 100 * change / mean( productivity_sim ) ;
    close all; 
    
    h = figure('position', [0, 0, 400, 260]);
    plot( lag_years,   smooth(ratio*coef{2}), 'LineStyle', '--', 'Color', [0.9100    0.4100    0.1700],  'LineWidth',1 ); hold on;  % Crisis
    plot( lag_years,   smooth(ratio*coef{3}), 'k-.' , 'LineWidth',1 );      % Recession
    plot(lag_years, 0*lag_years,':' );
    if( type=='Post' )  location = 'SouthWest';  else  location = 'NorthWest';      end
    legend(    [  "Crisis", "Recession"  ] , 'Location', location  ); legend boxoff ;
%     if( N_output>1 )
%         saveTightFigure( h,  [ 'ImpulseSummary_',  cell2mat(type)  ,'Crisis.pdf' ]  ); 
%     end
end


if( N_output>1 )
%% Part 2: Recovery Dynamics in Productivity and Credit Spreads
post_index = [ -5/dt: 10/dt ];  % -3 years to +10 years
GDP_recovery_crisis = zeros( numel(post_index), 1  ) ;         Spread_recovery_crisis = zeros( numel(post_index), 1  ) ; 
GDP_recovery_recession = zeros( numel(post_index), 1  ) ;   Spread_recovery_recession = zeros( numel(post_index), 1  ) ; 
InBuffer = (1:N_sim > 15/dt)  &  1:N_sim < (N_sim - 15/dt);
InBuffer = InBuffer';
for(  i = post_index  )
    GDP_recovery_crisis( i-min(post_index)+1 ) = mean(  productivity_sim( find(CrisisStart & InBuffer ) + i )    ) ;
    Spread_recovery_crisis( i-min(post_index)+1 ) = mean(  credit_spread_sim( find(CrisisStart& InBuffer) + i )    ) ;
    GDP_recovery_recession( i-min(post_index)+1 ) = mean(  productivity_sim( find(RecessionStart& InBuffer) + i )    ) ;
    Spread_recovery_recession( i-min(post_index)+1 ) = mean(  credit_spread_sim( find(RecessionStart& InBuffer) + i )    ) ;
end
types = [ "Productivity", "Credit Spread"];
for( k = 1:numel(types) )
    type=types(k);
    if( type== "Productivity" )
        avg= mean(productivity_sim);    crisis_quantity = GDP_recovery_crisis;    recession_quantity = GDP_recovery_recession;
    else
        avg= mean(credit_spread_sim);    crisis_quantity = Spread_recovery_crisis;    recession_quantity = Spread_recovery_recession;
    end
    h = figure('position', [0, 0, 400, 300]);
    plot(  post_index/dt,  100*(crisis_quantity - crisis_quantity( post_index==-1 ) ) / avg ,   post_index/dt, 100*( recession_quantity - recession_quantity(post_index==-1) ) / avg, '--'  )
    xlabel(  'Post Crisis/Recession Year' );   ylabel( [ cell2mat(type) , 'Percentage Change' ] );  legend( [ "Crisis", "Recession"  ]  );   
%     saveTightFigure( h,  [  cell2mat(type),  ' Recovery_model simulation.pdf'  ] )
end


%% Illustration: The Surprise Component of Bank-Run. Larger Surprise, larger disruption
T_recovery = 10; % Simulations for 10 years
dt = 1/52; % The unit of simulation is 1 week. 
N_recovery = T_recovery / dt; 
Years = (1:N_recovery)*dt;
w_mean = mean( w_sim );

w1_benchmark = mean(w_sim(RecessionStart)) ;    % Starting from a recession
w1 = w1_benchmark - ppval( kappaw_fitted, w1_benchmark ); 

w2_benchmark = mean(w_sim(CrisisStart));   % Starting from a crisis
%lambda2 = lambda_mean + kappa_lambda_fun( lambda_mean );
w2 = w2_benchmark - ppval( kappaw_fitted, w2_benchmark ); 

productivity_path1 = zeros( N_recovery, 1 ); 
productivity_path1_benchmark = zeros( N_recovery, 1 ); 
productivity_path2 = zeros( N_recovery, 1 ); 
productivity_path2_benchmark = zeros( N_recovery, 1 ); 
dN=0;  dB=0;    % No shock after the initial shock

for(i=1:N_recovery)   
    productivity_path1(i) = ppval( productivity_fitted, w1 );
    productivity_path1_benchmark(i) = ppval( productivity_fitted, w1_benchmark);   % Note: the evolution of lambda is the same for both benchmark and the shock case   
    productivity_path2(i) = ppval( productivity_fitted, w2 );
    productivity_path2_benchmark(i) = ppval( productivity_fitted, w2_benchmark );   % Note: the evolution of lambda is the same for both benchmark and the shock case   
    
    muw1 = ppval( muw_fitted, w1 );  
    muw1_benchmark =  ppval( muw_fitted, w1_benchmark );  
    muw2 = ppval( muw_fitted, w2 );  
    muw2_benchmark = ppval( muw_fitted, w2_benchmark );  
    
    w1 = w1 + muw1*dt  ;
    w1_benchmark = w1_benchmark + muw1_benchmark*dt  ;
    w2 = w2 + muw2*dt  ;
    w2_benchmark = w2_benchmark + muw2_benchmark*dt  ;
end

close all; h = figure('position', [0, 0, 800, 300]);
index = 1:12:N_recovery;
subplot(1,2,1); plot( Years(index),  productivity_path1_benchmark(index), Years(index),  productivity_path1(index),   'bo',   Years(index),  productivity_path2_benchmark(index), '--', Years(index),   productivity_path2(index) ,  'o', 'MarkerSize', 2 ); 
legend( 'Recession',  'Bank-run After Recession', 'Crisis', 'Bank-run After Crisis', 'Location', 'southeast' );  title( 'Recovery Dynamics' )
subplot(1,2,2); plot( Years, smooth( log( productivity_path1./productivity_path1_benchmark )) ,   Years,  smooth( log( productivity_path2./productivity_path2_benchmark )) , '--' ); 
legend( 'Bank-Run After Recession', 'Bank-Run After Crisis', 'Location', 'southeast' );   title( 'Relative Recovery' );
% saveTightFigure(  h,   'BankRun_Surprise.pdf' )

simulation_results.w_limit = w_limit;
simulation_results.productivity_limit = productivity_limit;

end






